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Abstract 

Monte Carlo simulations of mono- and polydisperse two-dimensional crystals are reported. 
The particles in the studied system, interacting through hard-core repulsive Yukawa potential, 
form a solid phase of hexagonal lattice. The elastic properties of crystalline Yukawa systems are 
determined in the NpT ensemble with variable shape of the periodic box. Effects of the Debye 
screening length contact value of the potential (e), and the size polydispersity of particles 

on elastic properties of the system are studied. The simulations show that the polydispersity 
of particles strongly influences the elastic properties of the studied system, especially on the 
shear modulus. It is also found that the elastic moduli increase with density and their growth 
rate depends on the screening length. Shorter screening length leads to faster increase of elastic 
moduli with density and decrease of the Poisson’s ratio. In contrast to its three-dimensional 
version, the studied system is non-auxetic, i.e. shows positive Poisson’s ratio. 


1 Introduction 

Colloidal suspensions are systems consisting of small particles of sizes from tens to hundreds nanome¬ 
ters dispersed in a solvent [l|. The particles can form colloidal crystals not only in three dimensions 
(3D)i but also in two dimensions (2D) where they create a hexagonal lattice j^. When particles are 
endowed with electrical charge and placed in a dispersing mediumUhey exhibit exponential decay of 
electric potential due to screening effect from ions in the solvent (J]. This decay is described by the 
screening length k~^ (known also as the Debye length) which depends on temperature, permittivity, 
and concentration of ions in the solvent. Interactions between particles in such systems are decribed 
by Hard-Core Repulsive Yukawa Potential (HCRYP) [1, @|, both in two Q and three d, @| dimen¬ 
sions. The potential consists of two parts: the hard core (reflecting the influence of finite sizes of 
particles) and the repulsive part describing spherically symmetric Coulomb screening effect. HCRYP 
systems can behave (d^ending on k~^) as hard particle systems (perfect screening) or as Wigner 
crystal (no screening) [9|. 

More than ten years ago the liquid-solid phase diagrams for two-dimensional (2D) j3| and three- 
dimensional (3D) 3 3, m H HCRYP systems have been studied. The research showed that in 
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Figure 1: Size dependence of the elastic moduli B (a), ^ (b), Poisson’s ratio (c), and Young’s 
modulus (d) for the system consisting oi N = 56, 224, 504, 896, 1400 particles, at p* = 30, /3e = 20 
and KC = 10. The lines represent linear fits to the obtained data. 


2D HCRYP particles can crystallize into a hexagonal lattice, whereas the in 3D they form f.c.c. 
or b.c.c. structures (depending on external conditions). Various thermodynamic properties and 
physical phenomena have been studied for the Yukawa potential, both in 2D 0, H, [i3-[l3l and 3D 0- 


Ifl) 111 18-^. Recently the elastic properties of three dimensional monodisperse Yukawa system 
were studied and it was shown that the system is partially auxetic 25l |. 

The term auxetic ( 2 ^ refers to the class of materials and models that exhibit unusual elastic 
properties. They expand (shrink) transversally when stretched (compressed) longitudinally. Their 
Poisson’s ratio [29| (in contrast to common materials) is negative. Due to their exceptional properties, 
auxetic materials and models gather ever growing interest in recent years j^. 

It is worth noting that usually colloidal systems are not composed of identical particles. Thus 
it is meaningful to consider the enhancement of the model by introduction of the particle size dis¬ 
persion 2l|, [ 2 ^. The study shows that the dispersion in sizes of particles can significantly impact 


the thermodynamic |3ll . l32l | (especially elastic |33l | ) properties of the system. The study of poly- 


disperse hard disks system not only indicated a strong dependence of elastic constants on particles 
size pojydispersity, but also revealed surprising behaviour of the Poisson’s ratio in the close-packing 


limit 3J|. In this context, the subject of this work are both mono- and poly disperse systems. 


This work encompasses the results of study performed on two dimensional hexagonal systems 
composed of Yukawa particles. The elastic properties of mono- and polydisperse two-dimensional 
crystals of hard-core repulsive Yukawa particles are described. Beside studies of the influence of 
screening length, the temperature and the particles size polydispersity on elastic properties, the 
authors posed additional question: whether this system will also exhibit auxetic properties. It is 
worth to add that the elastic properties of such a system have never been investigated. 

The structure of this manuscript is the following. In section [2] basic formulas and definitions of 
physical quantities necessary to describe the elasticity of the studied system, as well as the description 
of models and methods are provided. In section |3] the simulation results and their discussion are 
presented. The last section (jT]) provides the summary of the work. 


2 Preliminaries 

2.1 Elastic properties 

of the system under non-zero external pressure p which is subjected to reversible deformation can 
be described by expanding the free enthalpy (Gibbs free energy) with respect to the Lagrange strain 
tensor components. The change of the Gibbs free energy corresponding the elastic deformation of 
hexagonal lattice under external pressure p expressed with respect to the macroscopic elastic moduli. 
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Figure 2: Thermodynamic characteristics of the monodisprese hard-core repulsive Yukawa system 
at KU = 10. The pressure p*(a); bulk modulus i?*(b); shear modulus /U*(c) and the Poisson’s ratio 
i^(d) are shown as functions of density for three selected values of temperatures. The lines are drawn 
to guide the eye. 



Figure 3: Thermodynamic characteristic of the monodisperse hard-core repulsive Yukawa system at 
/3e = 20 as functions of density: (a) pressure p*, (b) the bulk modulus B*, (c) the shear modulus /r*, 
and (d) the Poisson’s ratio. The lines are drawn to guide the eye. 


is of the form: 


AG/Yp = 


+ 


-B{exx + £yyf 

—+ {Exx ~ £yy) ] ) 


( 1 ) 


where e is the strain tensor; Vp is the reference volume (at equilibrium) of the system under pressure 
p; B is the bulk modulus, and u is the shear modulus (s^. It follows from eq. ([T]) that hexagonal 
lattice is elastically isotropic (29|. 


2.2 The Poisson’s ratio 


is the negative ratio of the relative changes of transverse to longitudinal dimensions when only the 
longitudinal stress tensor component is infinitesimally changed. For hexagonal structure in two 
dimensions under pressure p one can express this quantity in the form 3^ : 

B — pL 


V = 


B p 

It is worth adding that for 2D systems the Poisson’s ratio can change in the range of — I < < I 


( 2 ) 


2.3 Particle size dispersion 

was introduced randomly according to the normal distribution with the first two central moments 
equal to a and <5, respectively. The standard deviation of the particle size distribution divided by 
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the mean diameter of the particle, called the polydispersity parameter [32|: 


1 


^ = -r\\l (o-^) - > 

(cr) V 


( 3 ) 


was used to control the degree of disorder in the system. It should be noted that, since the system 
consists of finite number of particles, it is important to control the values of {a) and 5 as well 
as possible. The simplest way to meet this condition is to repeatedly generate different sets of 
random numbers until the one is found whose properties fall, within a specified accuracy, near the 
desired values (3fi| . In this work, however, a different approach was used. After generating a set 
of n Gaussian pseudo random numbers (using method described in |37|), the corresponding average 
and standard deviation are used to convert this set to the standard normal distribution 
Z = {X — . Next one can convert obtained values such that they exactly match any desired 

distribution X' = (Z — a)/6, with the machine accuracy. This approach ensures the compatibility 
of distributions for different samples with very high accuracy. 


2.4 The model description 


starts with introducing reduced thermodynamic quantities used in the simulations: 


reduced value 

symbol 

therm, equiv. 

inverse temperature 

/3e 

elksT 

number density 

P* 

pa^ 

pressure 

P* 

pPa'^ 

bulk modulus 

B* 

Ba^/ksT 

shear modulus 

P* 

pa’^/ksT 


where T is the temperature and ks is the Boltzman constant. The studied system consists of N 
particles interacting via HCRYP. For the monodispersed systems, the potential of interaction is of 
the form: 




oo, 

Q, exp[-K(r-o-)] 

rju 


r < a, 
r > a, 


(4) 


where r is the distance between the centers of particles, a is the diameter of the hard core, e is 
the value of pair interaction of particles at contact (i.e. at r = cr) and is the Debye screening 
length. In the case of systems with non-zero size polydispersity of particles a generalization of the 
equation (j4]) proposed by Colombo and Dijkstra [2^ was used: 


/3u{r) 


oo, 

/3e^exp[-K(r-ay)], 


r < aij, 
r > aij, 


(5) 


2.5 Computer simulations 

have been performed in the isothermal-isobaric ensemble (NpT) by Monte Carlo method. In order to 
calculate the elastic moduli, the approach based on the averaging of strain fluctuations proposed by 
Parrinello-Rahman and further developed by Ray and Rahman 39, |^, was used. The details 


of the method were described in ref. |4l| . It is worth notin g th at in the case of a potential with 


long-range interactions the minimum image method (MIM) 


can be used. Our studies showed 
that to obtain accurate results for following parameters (k < 4, /3 > 80, rcut < 2.5) of studied 
potential the MIM or another method, which takes into account the long-range interactions, must 
be used. 


4 












The typical studied systems were composed of = 224 particles under periodic boundary con¬ 
ditions. Samples for simulation were chosen such that the edge ratio was as close to 1 as possible. 
The particles were placed in the 14a x 8^/3a box, where a is the distance between the nearest 
neighbors in crystal. The term nearest neighbor refers to particles that share a common edge of 
their Dirichlet polygons. Results from simulations were averaged over 10^ Monte Carlo cycles af¬ 
ter prior equilibration of 10® cycles. To check the dependence of the elastic moduli on the size of 
studied systems, the system of = 56, 224, 504, 896, 1400 particles have been studied, for which 
5 • 10®, 10^, 3 • 10^, 6 • 10^, 10® MC cycles were performed. The interaction potential implemented 
for simulations of polydisperse systems is described by the equation ([5]). For monodisperse systems 
the latter reduces to the form of eq. (HD- Elastic moduli for each phase point were determined by 
averaging over 10 independent simulations, each of which was performed for a different structure 
generated from a proper distribution. 

3 Results and discussion 
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Figure 4: Equation of state of the hard-core repulsive Yukawa system at /3e = 20 for kct = 5 (a), 
Ka = 10 (b), and na = 20 (c). The lines are drawn to guide the eye. 



0.6 0.8 1 
P* 


1.4 

1.2 




0.8 


0.6 


X = fi* B* Ka = 10 
5 = 0.02 • O 
S = 0.04 A A 
5 = 0.06 ♦ O 
5 = 0.08 ▼ V 




0.6 


0.8 

P* 


'O 


1.4 - 

1.2 - 

1 

0.8 

0.6 - 


X = iT B' 
5 = 0.02 • O 
5 = 0.04 A A 

4 = 0.06 ♦ * 

5 = 0.08 T V 


Ku — 20 


0.8 


Figure 5: Relative change of elastic moduli due to the non-zero size polydispersity for kcj = 5 (a), 
Ka = 10 (b) and Ka = 20 (c). X stands for B or The lines are drawn to guide the eye. 


Based on the equation of state in Ref. [3|, values of the potential’s parameters corresponding to 
solid phases were chosen to allow for comparison of simulation results with experimental data, e.g. 
/3e ~ 20 corresponds to colloidal systems in a low dielectric solvent at room temperature |22l |. The 
discussion of elastic properties of given systems that follows, is presented in the form of the three 
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Figure 6: Relative change of (a) the bulk modulus and (b) the shear modulus due to the non-zero 
size polydispersity, studied at (3e = 20, for two different values of pressure. The lines are drawn to 
guide the eye. 


cases: the constant screening length (ko") ^ (sec. EU; the constant temperature (/3e) ^ Isec. 13.2)1 : 
and the influence of particle size polydispersity (sec. ESI) on the mentioned above properties. 

In Figure [T] are shown dependencies of the elastic moduli and Poisson’s ratio with respect to the 
size of the system. It may be noted that the simulations performed for N = 224 yield results that 
differ only a few percent from the values extrapolated to the thermodynamic limit (N —)• oo). One 
can add that the agreement is quite good for all the studied system sizes. In view of the above, all 
further discussion applies to systems containing N = 224 particles. 

3.1 The case of constant screening length 

is illustrated in Figure El The elastic moduli as functions of density for («:cr)“^ = 10“^ are presented 
there. One can observe a monotonic increase of both elastic moduli and pressure with increasing 
density. A decrease of the temperature results in an increase of pressure and elastic moduli at entire 
range of studied densities. Furthermore, it causes a shift of the entire isotherm describing the lower 
temperature in the direction of high densities. 

The density dependence of Poisson’s ratio is presented in Fig. Eh- h is immediately apparent 
that for selected n values of Poisson’s ratio for different temperatures are very close one to another 
at high densities. This suggests that the influence of particles’ hard cores becomes dominant at 
these densities. At low densities an increase of temperature causes an increase of the Poisson’s ratio, 
however the dependence on temperature is rather weak. 

3.2 The case of constant temperatnre 

was investigated for = 20“^. Figure El illustrates the density dependence of pressure (a) 

as well as the bulk modulus (b), the shear modulus (c) and the Poisson’s ratio (d). Similarly to 
previous case, an increase of the pressure and the elastic moduli with the density can be observed. 
However, the growth rate of the elastic moduli is different and it clearly depends on the screening 
length. An increase of the elastic moduli with increasing density is slower for long-range potential 
((«:cj)“^ = 5“^) than for interactions with shorter screening length ((«:cr)“^ = 20“^). 

Figure Eh illustrates the dependence of Poisson’s ratio on density for systems with different 
screening lengths. It is noticeable that for a given k the Poisson’s ratio weakly depends on the density 
of the system. However, the screening length itself essentially influences the value of Poisson’s ratio. 
This leads to conclusion that one can control the Poisson’s ratio by modifying the Debye screening 
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Figure 7: Density dependence of the Poisson’s ratio of the polydisperse hard-core repulsive Yukawa 
system at /3e = 20, for three values of the screening length kcj = 5 (a), Ka = 10 (b), and Ka = 20 (c). 


length (e.g. by changing the properties of the electrolyte). 

3.3 Influence of size polydispersity 

of the particles’ hard cores on elastic properties is discussed in this section. One can start from the 
examination of the equation of state presented in the Figure U) As in the previous case, presented 
data correspond to systems at reduced temperature (/3e)“^ = 20~^ for three values of screening 
length 5“^, 10“^, and 20“^. The choice of these values is dictated by two factors. Firstly, 

studies of Poisson’s ratio in case of monodisperse system (see sec. EU showed that its temperature 
dependence is rather week while the screening length has a significant impact on its value. Secondly, 
the present study is mainly focused on looking for mechanisms affecting Poisson’s ratio. Thus the 
influence of the screening length on the elastic properties and Poisson’s ratio in systems with non-zero 
polydispersity of particles sizes was considered at a selected (constant) temperature. Systems with 
four different values of polydispersity parameter 6 have been investigated. Its influence manifests 
as the change of number density (i.e. the ratio of the number of particles per 2D ’volume’) for a 
given pressure (see Fig. HD. An increase of the size polydispersity causes a decrease of the number 
density at a given pressure. The decrease, however, depends on the screening length. For long-range 
interactions (Fig. 0^) the impact of size dispersion is significantly smaller and noticeable only at high 
pressures. Once the screening length is getting smaller (increasing k), rising 5 enforces the system 
to expand even at lower pressures. 

Figure [5] illustrates changes of elastic moduli caused by the size polydispersity relative to monodis¬ 
perse system for different screening lengths. The relations are presented as functions of density. One 
can notice that the influence of the size polydispersity is very different for both the moduli. In the 
case of the bulk modulus, a signihcant increase can be observed only for high densities and long 
interactions (Fig[5^). For lower densities the influence of the size polydispersity on the bulk modulus 
is negligible. Furthermore, when the effect of particles’ hard cores becomes dominant (increasing 
k), one can observe only minor impact of the non-zero polydispersity on B in the whole range of 
studied densities (Fig. [5t;). This effect is further illustrated in the Figure [6^, where the relative bulk 
modulus with respect the polydispersity parameter is presented for three screening lengths at high 
and low pressures. It can be seen there that only for long-range interactions and high pressure one 
can observe a significant influence of polydispersity on the bulk modulus, and this impact quickly 
decays with the decrease of the screening length. 

The effect of the size polydispersity on the shear modulus /r is stronger than on the bulk mod¬ 
ulus. In Figure [5^ one can see that an increase of 6 decreases the shear modulus. From the other 
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hand shortening of the screening length increases an influence of the polydispersity on the shear 
modulus (Fig. [5] and [6]). Thus, one can conclude that for a system with given size polydispersity 
of particles one can additionally control its (the system) susceptibility to the shear deformations by 
changing the Debye screening length (e.g. by modifying the properties of the electrolyte). Unlike 
the bulk modulus, pressure has a little impact on a decrease of /r with the increasing parameter of 
polydispersity (Fig. [hja). Only in the case of shortest screening lengths (kct = 20) the difference 
between the considered properties of systems under low and high pressures is noticeable. For smaller 
values of k only an increase of the polydispersity parameter is responsible for a decrease of the shear 
modulus. It is worth noting, that most significant decrease of the shear modulus is observed for the 
shortest screening length. 

Figure [7] shows the density dependence of the Poisson’s ratio for different values of polydispersity 
parameter and various screening lengths. It can be seen that for long-range interactions the size 
polydispersity of particles affects weakly the Poisson’s ratio and it weakly depends on the density. 
With shortening the screening length, the impact of polydispersity becomes more noticeable, getting 
very significant role for the short-range interactions. 

4 Summary and conclusions 

The influence of the temperature, screening lengths and size polydispersity on elastic properties of 
the two-dimensional hexagonal structure of hard-core, repulsive Yukawa particles was investigated by 
Monte Carlo simulations in the NpT ensemble with variable shape of the periodic box. The obtained 
results show that influence of the screening length on the elastic properties of the studied systems 
is significant. It was also shown that size polydispersity impacts the elastic properties, especially 
the shear modulus, which decreases with decreasing the size polydispersity parameter. However, an 
increase of particle size dispersion is reflected in the increase of Poisson’s ratio. It was shown that 
the 2D Yukawa model is non-auxetic but the value of the Poisson’s ratio decreases with decreasing 
screening length. This mechanism can be useful in designing systems with low Poisson’s ratio. 
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